library(bugcount)
data_dir <- file.path(system.file(package = "bugcount"),"extdata")
Whitefly analysis.
Read WF count data.
WF counts per picture. Each picture must have a unique combination of
experiment + replicate + plant + pic
experiment description consists of:
exp_year + exp_cross + exp_propagation + exp_substrate
wf_file <- file.path(data_dir, "wf_consolidated.tab")
wf <- read_wf(file = wf_file, sep = "\t",header=TRUE)
Counts per Plant
wf_count <- plant_count(wf)
summary(wf_count$nymphs)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0 555 1859 2506 3460 25317
Internal check analysis.
wf_check <- wf_count[grepl("internal_check", wf_count$group) &
wf_count$clone != "Secundina",]
wf_check <- remove_levels(wf_check)
checks_count_fit <- count_fit(wf_check)
plot_count_fit(checks_count_fit, main = "Nymphs")

Assume that nymphs per plant are distributed as negative binomial
by_experiment <- fit_nb_glm(nymphs ~ experiment, wf_check)
plot_fit_nb_glm(wf_check,by_experiment, type = "violin", cld = FALSE)

add infestation regime as factor high,low: inferred from nyphs ~ experiment GLM “see below”
wf_check$infestation <- "low"
wf_check[wf_check$exp_year > 2016,]$infestation <- "high"
wf_check$infestation <- factor(wf_check$infestation, levels =c("low","high"))
Correlation between experiments
What indicator should we use as resistance phenotype? Nymph counts are not normal Error is not distributed normally. Nymph counts do not meet linear model assumptions (ANOVA, MapQTL). Best fit is a negative Binomial Distribution
cor_title <- paste("Geometric Means R2","Log Scale", sep = "\n")
check_by_exp <- wf_wide(clone ~ experiment, wf_check,
fun = function(x) log10(geometric_mean(x)) )
plot_check_cor(check_by_exp, main = paste("Checks", cor_title))

# nymphs ~ clone GLM posthoc and common letter difference
fit <- fit_nb_glm(nymphs ~ clone, wf_check)
plot_fit_nb_glm(wf_check,fit,type = "violin", xmax = 10000)

fit <- MASS::glm.nb(nymphs ~ infestation * clone, data = wf_check)
posthoc <- lsmeans::cld(lsmeans::lsmeans(fit, ~ infestation * clone,
adjust = "tuckey"), type = "response")
# Plot post hoc
### Order the levels for printing
clone_order <- wf_aggregate(wf_check, by_stat = "mean")["clone"]
posthoc$clone <- factor(posthoc$clone, levels = clone_order$clone)
posthoc$infestation = factor(posthoc$infestation,
levels=c("low", "high"))
### Remove spaces in .group
posthoc$.group=gsub(" ", "", posthoc$.group)
# quartz()
plot_gg_infestationXclone(posthoc)

Whole population.
Probability density model fit to actual distribution
wf_count_fit <- count_fit(wf_count)
plot_count_fit(wf_count_fit, main = "Nymphs")

Reproducibility
Do insect counts change per experiment?
nymphs ~ experiment GLM, posthoc and common letter difference
by_experiment <- fit_nb_glm(nymphs ~ experiment, wf_count)
plot_fit_nb_glm(wf_count,by_experiment, type = "violin", cld = FALSE)
add infestation regime as factor high,low: inferred from nyphs ~ experiment GLM “see below”
wf_count$infestation <- "low"
wf_count[wf_count$exp_year > 2016,]$infestation <- "high"
wf_count$infestation <- factor(wf_count$infestation, levels =c("low","high"))
Correlation between experiments
cor_title <- paste("Geometric Means Correlation","Log Scale", sep = "\n")
for (cross in levels(wf_count$exp_cross)) {
wf_by_exp <- wf_wide(clone ~ experiment,
wf_count[wf_count$exp_cross == cross,])
plot_wide_cor(wf_by_exp, main = paste(cross,cor_title))
}


Nymph counts by Cross

LS0tCnRpdGxlOiAiQnVnIENvdW50IEFuYWx5c2lzIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgIGhpZ2hsaWdodDogdGFuZ28KICAgIG51bWJlcl9zZWN0aW9uczogeWVzCiAgICB0aGVtZTogc3BhY2VsYWIKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwotLS0KCmBgYHtyfQpsaWJyYXJ5KGJ1Z2NvdW50KQpkYXRhX2RpciA8LSBmaWxlLnBhdGgoc3lzdGVtLmZpbGUocGFja2FnZSA9ICJidWdjb3VudCIpLCJleHRkYXRhIikKCmBgYAoKIyBXaGl0ZWZseSBhbmFseXNpcy4KIyMgUmVhZCBXRiBjb3VudCBkYXRhLiAKV0YgY291bnRzIHBlciBwaWN0dXJlLiBFYWNoIHBpY3R1cmUgbXVzdCBoYXZlIGEgdW5pcXVlIGNvbWJpbmF0aW9uIG9mICAKZXhwZXJpbWVudCArIHJlcGxpY2F0ZSArIHBsYW50ICsgcGljICAKZXhwZXJpbWVudCBkZXNjcmlwdGlvbiBjb25zaXN0cyBvZjogIApleHBfeWVhciArIGV4cF9jcm9zcyArIGV4cF9wcm9wYWdhdGlvbiArIGV4cF9zdWJzdHJhdGUgIApgYGB7cn0Kd2ZfZmlsZSA8LSBmaWxlLnBhdGgoZGF0YV9kaXIsICJ3Zl9jb25zb2xpZGF0ZWQudGFiIikKd2YgPC0gcmVhZF93ZihmaWxlID0gd2ZfZmlsZSwgc2VwID0gIlx0IixoZWFkZXI9VFJVRSkKYGBgCiMjIENvdW50cyBwZXIgUGxhbnQgIAoKYGBge3J9Cgp3Zl9jb3VudCA8LSBwbGFudF9jb3VudCh3ZikKc3VtbWFyeSh3Zl9jb3VudCRueW1waHMpCgpgYGAKCiMjIEludGVybmFsIGNoZWNrIGFuYWx5c2lzLgoKYGBge3IsIHdhcm5pbmc9RkFMU0V9CndmX2NoZWNrIDwtIHdmX2NvdW50W2dyZXBsKCJpbnRlcm5hbF9jaGVjayIsIHdmX2NvdW50JGdyb3VwKSAmCiAgICAgICAgICAgICAgICAgICAgICAgd2ZfY291bnQkY2xvbmUgIT0gIlNlY3VuZGluYSIsXQoKd2ZfY2hlY2sgPC0gcmVtb3ZlX2xldmVscyh3Zl9jaGVjaykKCmNoZWNrc19jb3VudF9maXQgPC0gIGNvdW50X2ZpdCh3Zl9jaGVjaykKCnBsb3RfY291bnRfZml0KGNoZWNrc19jb3VudF9maXQsIG1haW4gPSAiTnltcGhzIikKYGBgCgpBc3N1bWUgdGhhdCBueW1waHMgcGVyIHBsYW50IGFyZSBkaXN0cmlidXRlZCBhcyBuZWdhdGl2ZSBiaW5vbWlhbAoKYGBge3J9CmJ5X2V4cGVyaW1lbnQgPC0gZml0X25iX2dsbShueW1waHMgfiBleHBlcmltZW50LCB3Zl9jaGVjaykKCnBsb3RfZml0X25iX2dsbSh3Zl9jaGVjayxieV9leHBlcmltZW50LCB0eXBlID0gInZpb2xpbiIsIGNsZCA9IEZBTFNFKQpgYGAKYWRkIGluZmVzdGF0aW9uIHJlZ2ltZSBhcyBmYWN0b3IgCmhpZ2gsbG93OiBpbmZlcnJlZCBmcm9tIG55cGhzIH4gZXhwZXJpbWVudCBHTE0gInNlZSBiZWxvdyIKCmBgYHtyfQp3Zl9jaGVjayRpbmZlc3RhdGlvbiA8LSAibG93Igp3Zl9jaGVja1t3Zl9jaGVjayRleHBfeWVhciA+IDIwMTYsXSRpbmZlc3RhdGlvbiA8LSAiaGlnaCIKd2ZfY2hlY2skaW5mZXN0YXRpb24gPC0gZmFjdG9yKHdmX2NoZWNrJGluZmVzdGF0aW9uLCBsZXZlbHMgPWMoImxvdyIsImhpZ2giKSkKYGBgCgoKIyMjIENvcnJlbGF0aW9uICBiZXR3ZWVuIGV4cGVyaW1lbnRzCldoYXQgaW5kaWNhdG9yIHNob3VsZCB3ZSB1c2UgYXMgcmVzaXN0YW5jZSBwaGVub3R5cGU/Ck55bXBoIGNvdW50cyBhcmUgbm90IG5vcm1hbCAKRXJyb3IgaXMgbm90IGRpc3RyaWJ1dGVkIG5vcm1hbGx5LgpOeW1waCBjb3VudHMgZG8gbm90IG1lZXQgbGluZWFyIG1vZGVsIGFzc3VtcHRpb25zIChBTk9WQSwgTWFwUVRMKS4KQmVzdCBmaXQgaXMgYSBuZWdhdGl2ZSBCaW5vbWlhbCBEaXN0cmlidXRpb24KCmBgYHtyLCBmaWcuYXNwPTF9CmNvcl90aXRsZSA8LSBwYXN0ZSgiR2VvbWV0cmljIE1lYW5zIFIyIiwiTG9nIFNjYWxlIiwgc2VwID0gIlxuIikKCmNoZWNrX2J5X2V4cCA8LSB3Zl93aWRlKGNsb25lIH4gZXhwZXJpbWVudCwgd2ZfY2hlY2ssIAogICAgICAgICAgICAgICAgICAgICAgICBmdW4gPSBmdW5jdGlvbih4KSBsb2cxMChnZW9tZXRyaWNfbWVhbih4KSkgKQpwbG90X2NoZWNrX2NvcihjaGVja19ieV9leHAsIG1haW4gPSBwYXN0ZSgiQ2hlY2tzIiwgY29yX3RpdGxlKSkKYGBgCgoKCmBgYHtyfQojIG55bXBocyB+IGNsb25lIEdMTSBwb3N0aG9jIGFuZCBjb21tb24gbGV0dGVyIGRpZmZlcmVuY2UKZml0IDwtIGZpdF9uYl9nbG0obnltcGhzIH4gY2xvbmUsIHdmX2NoZWNrKQpwbG90X2ZpdF9uYl9nbG0od2ZfY2hlY2ssZml0LHR5cGUgPSAidmlvbGluIiwgeG1heCA9IDEwMDAwKQpgYGAKCmBgYHtyLCBmaWcuYXNwPTAuNX0KZml0IDwtIE1BU1M6OmdsbS5uYihueW1waHMgfiAgaW5mZXN0YXRpb24gKiBjbG9uZSwgZGF0YSA9IHdmX2NoZWNrKQpwb3N0aG9jIDwtIGxzbWVhbnM6OmNsZChsc21lYW5zOjpsc21lYW5zKGZpdCwgfiBpbmZlc3RhdGlvbiAqIGNsb25lLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGFkanVzdCA9ICJ0dWNrZXkiKSwgdHlwZSA9ICJyZXNwb25zZSIpCiMgUGxvdCBwb3N0IGhvYwoKIyMjIE9yZGVyIHRoZSBsZXZlbHMgZm9yIHByaW50aW5nCmNsb25lX29yZGVyIDwtIHdmX2FnZ3JlZ2F0ZSh3Zl9jaGVjaywgYnlfc3RhdCA9ICJtZWFuIilbImNsb25lIl0KCnBvc3Rob2MkY2xvbmUgPC0gZmFjdG9yKHBvc3Rob2MkY2xvbmUsIGxldmVscyA9IGNsb25lX29yZGVyJGNsb25lKQoKcG9zdGhvYyRpbmZlc3RhdGlvbiA9IGZhY3Rvcihwb3N0aG9jJGluZmVzdGF0aW9uLAogICAgICAgICAgICAgICAgICAgbGV2ZWxzPWMoImxvdyIsICJoaWdoIikpCgojIyMgIFJlbW92ZSBzcGFjZXMgaW4gLmdyb3VwICAKCnBvc3Rob2MkLmdyb3VwPWdzdWIoIiAiLCAiIiwgcG9zdGhvYyQuZ3JvdXApCgojIHF1YXJ0eigpCnBsb3RfZ2dfaW5mZXN0YXRpb25YY2xvbmUocG9zdGhvYykKYGBgCgoKCiMjIFdob2xlIHBvcHVsYXRpb24uCgojIyMgUHJvYmFiaWxpdHkgZGVuc2l0eSBtb2RlbCBmaXQgdG8gYWN0dWFsIGRpc3RyaWJ1dGlvbgpgYGB7ciwgd2FybmluZz1GQUxTRX0Kd2ZfY291bnRfZml0IDwtICBjb3VudF9maXQod2ZfY291bnQpCgpwbG90X2NvdW50X2ZpdCh3Zl9jb3VudF9maXQsIG1haW4gPSAiTnltcGhzIikKYGBgCgojIyMgUmVwcm9kdWNpYmlsaXR5IAoKRG8gaW5zZWN0IGNvdW50cyBjaGFuZ2UgcGVyIGV4cGVyaW1lbnQ/CgpueW1waHMgfiBleHBlcmltZW50IEdMTSwgcG9zdGhvYyBhbmQgY29tbW9uIGxldHRlciBkaWZmZXJlbmNlIAoKYGBge3J9CmJ5X2V4cGVyaW1lbnQgPC0gZml0X25iX2dsbShueW1waHMgfiBleHBlcmltZW50LCB3Zl9jb3VudCkKCnBsb3RfZml0X25iX2dsbSh3Zl9jb3VudCxieV9leHBlcmltZW50LCB0eXBlID0gInZpb2xpbiIsIGNsZCA9IEZBTFNFKQpgYGAKYWRkIGluZmVzdGF0aW9uIHJlZ2ltZSBhcyBmYWN0b3IgCmhpZ2gsbG93OiBpbmZlcnJlZCBmcm9tIG55cGhzIH4gZXhwZXJpbWVudCBHTE0gInNlZSBiZWxvdyIKCmBgYHtyfQp3Zl9jb3VudCRpbmZlc3RhdGlvbiA8LSAibG93Igp3Zl9jb3VudFt3Zl9jb3VudCRleHBfeWVhciA+IDIwMTYsXSRpbmZlc3RhdGlvbiA8LSAiaGlnaCIKd2ZfY291bnQkaW5mZXN0YXRpb24gPC0gZmFjdG9yKHdmX2NvdW50JGluZmVzdGF0aW9uLCBsZXZlbHMgPWMoImxvdyIsImhpZ2giKSkKYGBgCiMjIyBDb3JyZWxhdGlvbiAgYmV0d2VlbiBleHBlcmltZW50cyAKCmBgYHtyLCBmaWcuYXNwPTF9CgoKZm9yIChjcm9zcyBpbiBsZXZlbHMod2ZfY291bnQkZXhwX2Nyb3NzKSkgewogIAogIHdmX2J5X2V4cCA8LSB3Zl93aWRlKGNsb25lIH4gZXhwZXJpbWVudCwKICAgICAgICAgICAgICAgICAgICAgICB3Zl9jb3VudFt3Zl9jb3VudCRleHBfY3Jvc3MgPT0gY3Jvc3MsXSkKICAKICBwbG90X3dpZGVfY29yKHdmX2J5X2V4cCwgbWFpbiA9IHBhc3RlKGNyb3NzLGNvcl90aXRsZSkpCiAgCn0KYGBgCgoKCiMjIyBOeW1waCBjb3VudHMgYnkgQ3Jvc3MgCmBgYHtyLCBmaWcuYXNwPTF9CmV4cF9sZXZlbHMgPC0gbGV2ZWxzKGFzLmZhY3Rvcih3Zl9jb3VudCRleHBlcmltZW50KSkKZXhwX2dycCA8LSBsaXN0KENNODk5NiA9IGxpc3QoMSwyLDQsNSksCiAgICAgICAgICAgICAgICBHTTg1ODYgPSBsaXN0KDMsNiw3KSkKcGxvdF9saXN0IDwtIGxpc3QoKQpmb3IgKGNyb3NzIGluIGxldmVscyh3Zl9jb3VudCRleHBfY3Jvc3MpKSB7CiAgZm9yIChpZHggaW4gMTpsZW5ndGgoZXhwX2dycFtbY3Jvc3NdXSkpIHsKICBleHBfYWxsb3dlZCA8LSBleHBfbGV2ZWxzW2V4cF9ncnBbW2Nyb3NzXV1bW2lkeF1dXQogIAogIHdmX2J5X2Nyb3NzIDwtIHdmX2NvdW50W3dmX2NvdW50JGV4cF9jcm9zcyA9PSBjcm9zcyAmCiAgICAgICAgICAgICAgICAgICAgICAgICAgICB3Zl9jb3VudCRleHBlcmltZW50ICVpbiUgZXhwX2FsbG93ZWQsXQogIGV4cF9jb3VudCA8LSBhZ2dyZWdhdGUoZXhwZXJpbWVudCB+IGNsb25lICx3Zl9ieV9jcm9zcywKICAgICAgICAgICAgICAgICAgICAgICAgIEZVTiA9IGZ1bmN0aW9uKHgpIGxlbmd0aCh1bmlxdWUoeCkpKQogIHNlbGVjdGVkX2Nsb25lcyA8LSB3aXRoKGV4cF9jb3VudCwgewogICAgY2xvbmVbZXhwZXJpbWVudCA9PSBsZW5ndGgoZXhwX2FsbG93ZWQpXQogIH0pCiAgCiAgd2ZfYnlfY3Jvc3MgPC0gd2ZfYnlfY3Jvc3NbIHdmX2J5X2Nyb3NzJGNsb25lICVpbiUgc2VsZWN0ZWRfY2xvbmVzLF0KICB3Zl9ieV9jcm9zcyA8LSAod2ZfYnlfY3Jvc3MpCiAgCiAgd2ZfYnlfY3Jvc3MkbnltcGhzIDwtIHdmX2J5X2Nyb3NzJG55bXBocysxCiAgIyBHTE0gKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqKioqICMKICAjIEFzc3VtaW5nIE5lZ2F0aXZlIEJpbm9taWFsIGRpc3RyaWJ1dGlvbgogICMKICAKICBmaXRfbmIgPC0gTUFTUzo6Z2xtLm5iKG55bXBocyB+ICBjbG9uZSwgZGF0YSA9IHdmX2J5X2Nyb3NzKQogIEFJQyhmaXRfbmIpCgogIHBvc3Rob2MgPC0gbHNtZWFuczo6Y2xkKAogIGxzbWVhbnM6OmxzbWVhbnMoCiAgICAgZml0X25iLCB+IGNsb25lLCBhZGp1c3QgPSAidHVja2V5IiksdHlwZSA9ICJyZXNwb25zZSIKICApCiAgCiAgIyBQbG90IHBvc3QgaG9jCiAgYnJlYWtzIDwtIHBvc3Rob2MkY2xvbmVbIWdyZXBsKCJHTXxDTSIsIHBvc3Rob2MkY2xvbmUsIHBlcmwgPSBUUlVFKV0KICAKICBwbG90X2xpc3RbW3Bhc3RlKGV4cF9hbGxvd2VkKV1dIDwtIHBsb3RfY3Jvc3NfbmJfZml0KHBvc3Rob2MsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBleHBfYWxsb3dlZCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJyZWFrcyA9IGJyZWFrcykKICAKICB9Cn0KZ2dwbG90Mjo6dGhlbWVfc2V0KGdncGxvdDI6OnRoZW1lX2dyYXkoYmFzZV9zaXplID0gMTApKQpvcmRlcmVkX2xpc3QgPC0gYyhwbG90X2xpc3RbMV0sbGlzdChlbXB0eSA9IE5VTEwpLHBsb3RfbGlzdFtjKDMsNywyLDUsNCw2KV0pCm5hbWVzKHBsb3RfbGlzdCkKbmFtZXMob3JkZXJlZF9saXN0KQpjb3dwbG90OjpwbG90X2dyaWQocGxvdGxpc3QgPSBvcmRlcmVkX2xpc3QsCiAgICAgICAgICBuY29sID0gMiwgbnJvdyA9IDQsIGFsaWduID0gInYiKQoKYGBgCg==